library(tidyverse)
library(lme4)
options(scipen=999)
df = read.csv(here::here('avg_tensor_by_roi_wide.csv'),
colClasses = c('subject' = 'factor',
'site' = 'factor',
'visit' = 'factor'))
outcomes = df %>%
select(where(is.numeric)) %>%
colnames
The purpose of the DTIMSCAMRAS study is to test for site effects on DTI metrics of different brain regions in the context of a traveling subjects design with a harmonized imaging protocol. Participants (n=11) with stable Multiple Sclerosis (MS) were scanned in 4 different locations: National Institute of Health (NIH), Brigham and Women’s Hospital (BWH), Johns Hopkins University (Hopkins), and The University of Pennsylvania (Penn). Two scans were collected per site visit, and participants were re-positioned between scans. Penn, BWH and NIH used Siemens scanners, while Hopkins used a Phillips scanner. The imaging modalities that were collected include T1s (MPRAGE), FLAIR, and Diffusion-weighted images (DWI).
To preprocess imaging data, qsiprep
was used for bias correction, skull-stripping and co-registration of the
T1s and DWIs, and specialized denoising and artifact correction for the
DWIs. Then, DTIFIT was used to compute the following scalar
maps:
Concurrently, several segmentation pipelines were used on the T1s (+ FLAIRs in the case of MIMOSA) with the purpose of defining particular regions of interest (ROIs):
To remove lesions appearing in GM or WM (or thalamus) from the calculation, mimosa masks were inverted and applied to the segmentation label maps in order to zero out the tissue masks at the location of each lesion. These adjusted segmentation results were then used to average the DTI scalar maps across all voxels belonging to each of the ROIs. The 36 resulting DTI metrics make up the set of outcome variables which are the focus of the site effects test. Below, the outcome variables are listed.
data.frame(OUTCOME = outcomes) %>%
separate(OUTCOME, c("SEGMENTATION", "ROI", "SCALAR"), remove = FALSE)
Note the general
Testing for site effects was performed using two approaches: a permutation-based approach and a deviance test of mixed models.
Additionally, as Hopkins data was acquired using a Philips scanner, which is in contrast to all other sites, the two testing approaches were performed again while excluding Hopkins data.
A permutation-based F-test was used to test for site effects. For each subject, the squared difference between all possible pairs of intra-site and inter-site measures were computed. For instance, take the first row of the subject data.frame below, which has BWH as the “reference site”: the squared differences between that row and the 6 non-BWH rows are computed for \(y\)——this process is repeated for all rows (and the results averaged) to arrive to the mean inter-site squared difference for this subject; to calculate the mean intra-site squared differences, the squared differences of the 4 intrasite pairs are averaged.
df %>%
filter(subject == '01001') %>%
select(!where(is.numeric), ATROPOS_GM_AD) %>%
rename(y = ATROPOS_GM_AD)
These differences were averaged across all subjects resulting in the Mean Squared difference Between Sites (\(\overline{MS_{B}}\)) and the Mean Squared difference Within Sites (\(\overline{MS_{W}}\)). The test statistic is composed of their ratio:
\[ F = \frac{\overline{MS_{B}}}{\overline{MS_{W}}} \]
For each metric, the observed test statistic was compared to a distribution of test statistics derived from 10000 permutations (i.e. a null distribution). The site labels were permuted within each subject. The proportion of null results higher than the observed test statistic served as a preliminary p-value \(p_0\). To prevent p-values of 0 (when the observed statistic was higher than all permuted results), the p-value was shifted using the following formula:
\[ p = \frac{10000*p_0 + 1}{10000 + 1}\]
As an additional test of site effects, mixed models were performed
corresponding to a crossed design using the lmer function
of the lme4 package. For each outcome variable, two
mixed-effects models were fitted: The formula for the base model
includes a random intercept for site as well as a random intercept for
subject. The formula for the extended model also includes these terms in
addition to a fixed effect term for site.
Base Model:
y ~ (1|subject) + (1|site)
Extended model:
y ~ site + (1|subject) + (1|site)
For each outcome, these two models were compared through a deviance
test using the anova function in R. P-values for these
tests are reported as a test of site effects.
Source code for this project can be found on GitHub.
Subjects 03001 and 03002 did not complete imaging at all 4 sites. Subject 03001 completed imaging at Penn and BWH only and 03002 completed imaging at BWH, Hopkins, and the NIH only.
Subject 02001 is missing DTI data from their NIH visit.
t_df <- df %>%
count(subject, site, .drop = FALSE) %>%
tidyr::pivot_wider(names_from = site, values_from = n) %>%
mutate(`Row Totals` = BWH + Hopkins + NIH + Penn)
t_df <- t_df %>%
summarize(across(where(is.integer), sum)) %>%
bind_rows(t_df, .) %>%
mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_df
Atropos segmentation failed for 6 images: 04001NIH01, 01003NIH01, 03002NIH02, 04003NIH01, 04001BWH02, and 03001BWH01.
knitr::include_graphics(here::here(c('misc/atropos_success.png', 'misc/atropos_fail.png')))
Left: Example of successful Atropos segmentation; Right: Example of a failed Atropos segmentation with one ROI missing
The following table excludes the failed atropos segmentations. Note subject 03001 only has 3 images after excluding failed atropos segmentations.
t_df <- df %>%
unite(sub, subject, site, visit, sep = '-') %>%
filter(!sub %in% c('04001-NIH-01', '01003-NIH-01', '03002-NIH-02', '04003-NIH-01', '04001-BWH-02', '03001-BWH-01')) %>%
separate(sub, c('subject', 'site', 'visit')) %>%
mutate_if(is.character, as.factor) %>%
count(subject, site, .drop = FALSE) %>%
tidyr::pivot_wider(names_from = site, values_from = n) %>%
mutate(`Row Totals` = BWH + Hopkins + NIH + Penn)
t_df <- t_df %>%
summarize(across(where(is.integer), sum)) %>%
bind_rows(t_df, .) %>%
mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_df
The failed atropos segmentations are included in the following visualizations but were excluded for the site effects tests.
TODO: compare harmonization with and without hopkins
Densities are colored by site, while the black line represents the overall density aggregated across sites. Distribution of the different sites tend to cluster together except for JLF Thalamus. Additionally, note the long tails caused by outliers in the ATROPOS distributions.
plot_densities <- function(var){
df %>%
ggplot(aes_string(x=var, color='site')) +
geom_density() +
stat_density(aes_string(x = var), geom = "line", color = "black") +
theme_bw() +
xlab(str_replace_all(var, "_", " "))
}
vars <-df %>%
select(ends_with('FA')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
Box plots show difference in medians across sites for most
variables.
Additionally, note the severe outliers in the ATROPOS metrics which
correspond to the failed segmentations. For FIRST THALAMUS, subject
01002 could be considered an outlier; their FA values are consistent
across sites and scans, however, suggesting this variation is meaningful
and that this subject should be retained.
plot_avg_tensor_values <- function(tensor){
df %>%
select(!is.numeric, ends_with(tensor)) %>%
pivot_longer(ends_with(tensor), names_to = 'seg', values_to = 'average') %>%
separate(seg, into = c('segmentation', 'roi', 'tensor')) %>%
unite(segmentation, segmentation, roi, sep = " ") %>%
mutate(segmentation = factor(segmentation,
levels = c("ATROPOS GM", "ATROPOS WM", "FAST GM", "FAST WM",
"JLF GM", "JLF WM", "JLF THALAMUS", "FIRST THALAMUS", "MIMOSA LESION"))) %>%
ggplot(aes(x=site, y=average)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(color = subject, shape=visit), alpha=0.5) +
#facet_grid(site~.) +
facet_grid(segmentation ~ .) +
coord_flip() +
ggtitle(sprintf("Mean %s in ROI by site", tensor)) +
theme_bw() +
# theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1)) +
xlab("") +
ylab("") +
theme(strip.text.y = element_text(angle = 0)) +
scale_x_discrete(expand = c(0.15, 0.15))
}
plot_avg_tensor_values('FA')
For Mean Diffusivity, distributions vary more widely across sites. In particular, MIMOSA LESION MD values are quite different for Hopkins data.
vars <-df %>%
select(ends_with('MD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
As before, medians appear to be different across sites.
plot_avg_tensor_values('MD')
Radial Diffusivity also shows some variations across sites (perhaps less than MD). JLF GM RD in NIH shows multiple peaks, and JLF WM RD shows skewed distributions for all sites. As before the metric for MIMOSA shows a different distribution for Hopkins data.
vars <-df %>%
select(ends_with('RD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
As before, medians appear to be differet across sites.
plot_avg_tensor_values('RD')
For Axial Diffusivity, Hopkins distributions show marked differences across most ROIs compared to other sites. The difference is quite pronounced for MIMOSA.
vars <-df %>%
select(ends_with('AD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
As before, medians show differences, particularly for Hopkins.
plot_avg_tensor_values('AD')
# load harmonized data
harmonized_df <- read.csv(here::here('results/avg_tensor_by_roi_wide_harmonized.csv'),
colClasses = c('subject' = 'character',
'site' = 'character',
'visit' = 'character'))
harmonized_no_hopkins_df <- read.csv(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_harmonized.csv'))
# load models
combat_model <- readRDS(here::here('results/avg_tensor_by_roi_wide_model.csv'))[c("gammahat", "delta2hat", "gammastarhat", "delta2starhat")]
combat_model_no_hopkins <- readRDS(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_model.csv'))[c("gammahat", "delta2hat", "gammastarhat", "delta2starhat")]
# model parameters
# what's the best way to plot combat estimates
Densities are colored by site, while the black line represents the overall density aggregated across sites. Distribution of the different sites tend to cluster together except for JLF Thalamus. Additionally, note the long tails caused by outliers in the ATROPOS distributions.
plot_densities <- function(var){
harmonized_df %>%
ggplot(aes_string(x=var, color='site')) +
geom_density() +
stat_density(aes_string(x = var), geom = "line", color = "black") +
theme_bw() +
xlab(str_replace_all(var, "_", " "))
}
vars <- harmonized_df %>%
select(ends_with('FA')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
Box plots show difference in medians across sites for most
variables.
Additionally, note the severe outliers in the ATROPOS metrics which
correspond to the failed segmentations. For FIRST THALAMUS, subject
01002 could be considered an outlier; their FA values are consistent
across sites and scans, however, suggesting this variation is meaningful
and that this subject should be retained.
plot_avg_tensor_values <- function(tensor){
harmonized_df %>%
select(!is.numeric, ends_with(tensor)) %>%
pivot_longer(ends_with(tensor), names_to = 'seg', values_to = 'average') %>%
separate(seg, into = c('segmentation', 'roi', 'tensor')) %>%
unite(segmentation, segmentation, roi, sep = " ") %>%
mutate(segmentation = factor(segmentation,
levels = c("ATROPOS GM", "ATROPOS WM", "FAST GM", "FAST WM",
"JLF GM", "JLF WM", "JLF THALAMUS", "FIRST THALAMUS", "MIMOSA LESION"))) %>%
ggplot(aes(x=site, y=average)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(color = subject, shape=visit), alpha=0.5) +
#facet_grid(site~.) +
facet_grid(segmentation ~ .) +
coord_flip() +
ggtitle(sprintf("Mean %s in ROI by site", tensor)) +
theme_bw() +
# theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1)) +
xlab("") +
ylab("") +
theme(strip.text.y = element_text(angle = 0)) +
scale_x_discrete(expand = c(0.15, 0.15))
}
plot_avg_tensor_values('FA')
For Mean Diffusivity, distributions vary more widely across sites. In particular, MIMOSA LESION MD values are quite different for Hopkins data.
vars <- harmonized_df %>%
select(ends_with('MD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
As before, medians appear to be different across sites.
plot_avg_tensor_values('MD')
Radial Diffusivity also shows some variations across sites (perhaps less than MD). JLF GM RD in NIH shows multiple peaks, and JLF WM RD shows skewed distributions for all sites. As before the metric for MIMOSA shows a different distribution for Hopkins data.
vars <-harmonized_df %>%
select(ends_with('RD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
As before, medians appear to be differet across sites.
plot_avg_tensor_values('RD')
For Axial Diffusivity, Hopkins distributions show marked differences across most ROIs compared to other sites. The difference is quite pronounced for MIMOSA.
vars <-harmonized_df %>%
select(ends_with('AD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)
for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}
As before, medians show differences, particularly for Hopkins.
plot_avg_tensor_values('AD')
# replace ATROPOS measures with NA for select images (segmentation failed)
fill_na_atropos <- function(df){
atropos_cols <- df %>%
select(contains('ATROPOS')) %>%
colnames()
df[(df$subject == '04001' & df$site == 'NIH' & df$visit == '01') |
(df$subject == '01003' & df$site == 'NIH' & df$visit == '01') |
(df$subject == '03002' & df$site == 'NIH' & df$visit == '02') |
(df$subject == '04003' & df$site == 'NIH' & df$visit == '01') |
(df$subject == '04001' & df$site == 'BWH' & df$visit == '02') |
(df$subject == '03001' & df$site == 'BWH' & df$visit == '01'), atropos_cols] <- NA
return(df)
}
df <- fill_na_atropos(df)
harmonized_df <- fill_na_atropos(harmonized_df)
The permutation-based F-test was carried out as previously described. The observed F-statistic for each of the metrics is plotted below as a black dot, while the distribution of permuted test statistics is shown as a white violin plot. All metrics showed significant site effects.
ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_null_F_df.rds")) %>%
mutate(outcome = map_chr(outcome, unique)) %>%
select(-rowname)
plot_F_dist <- function(ratio_stats, null_dist){
null_dist %>%
select(outcome, F) %>%
unnest() %>%
ggplot(aes(x=outcome, y=F)) +
geom_violin(draw_quantiles = c(0.95)) +
coord_flip() +
geom_point(aes(x=outcome, y=F),
data=tibble::rownames_to_column(ratio_stats))
}
plot_F_dist(ratio_stats, null_dist)
ratio_stats %>%
arrange(desc(F))
#' get p value of F by comparing to permutation distribution
test_ratio_stat = function(ratio.stats, null.dists){
p.vals = c()
for (col in colnames(ratio.stats)){
ratio.stat = pull(ratio.stats, col)
null.dist = pull(null.dists, col)
percentile = ecdf(null.dist)
p.val <- (10000* (1 - percentile(ratio.stat)) + 1)/(1+10000)
p.vals = c(p.vals, p.val)
}
names(p.vals) = colnames(ratio.stats)
return(p.vals)
}
#' make a data.frame with p values
make_perm_df <- function(ratio_stats, null_dist){
observed_F <- ratio_stats %>%
select(outcome, F) %>%
pivot_wider(names_from = 'outcome', values_from = 'F')
null_Fs <- null_dist %>%
select(outcome, F) %>%
unnest() %>%
pivot_wider(names_from = 'outcome', values_from = 'F') %>%
unnest() %>%
select(!!!colnames(observed_F)) # re-order columns to match observed
p_df <- data.frame(p_val = test_ratio_stat(observed_F, null_Fs))
perm_df <- cbind(
ratio_stats %>%
select(F),
p_df %>%
mutate(pval_fdr = p.adjust(p_val, method = 'fdr'))
)
return(perm_df)
}
perm_df <- make_perm_df(ratio_stats, null_dist)
plot_p_perm <- function(perm_df){
perm_df %>%
tibble::rownames_to_column('outcome') %>%
mutate(outcome = fct_reorder(outcome, pval_fdr, .desc = TRUE)) %>%
ggplot(aes(x=outcome, y=pval_fdr)) +
geom_point(size = 2) +
geom_hline(yintercept = 0.05, color = 'red') +
coord_flip() +
ggtitle("P-values of Permutation Test") +
ylim(0,1)
}
plot_p_perm(perm_df)
perm_df %>%
arrange(desc(p_val))
There were 35 out of 36 significant site effects according to permutation test.
Image acquisition at Hopkins was performed with a Philips scanner and distortion correction during on-scanner reconstruction. As the scanner manufacturer and reconstructions method differed than other sites, here we report the permutation results based on the non-Hopkins data only.
ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_null_F_df.rds")) %>%
mutate(outcome = map_chr(outcome, unique)) %>%
select(-rowname)
plot_F_dist(ratio_stats, null_dist)
ratio_stats %>%
arrange(desc(F))
perm_df <- make_perm_df(ratio_stats, null_dist)
plot_p_perm(perm_df)
perm_df %>%
arrange(desc(p_val))
ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_harmonized_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_harmonized_null_F_df.rds")) %>%
mutate(outcome = map_chr(outcome, unique)) %>%
select(-rowname)
plot_F_dist(ratio_stats, null_dist)
ratio_stats %>%
arrange(desc(F))
perm_df <- make_perm_df(ratio_stats, null_dist)
plot_p_perm(perm_df)
perm_df %>%
arrange(desc(p_val))
ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_harmonized_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_harmonized_null_F_df.rds")) %>%
mutate(outcome = map_chr(outcome, unique)) %>%
select(-rowname)
plot_F_dist(ratio_stats, null_dist)
ratio_stats %>%
arrange(desc(F))
perm_df <- make_perm_df(ratio_stats, null_dist)
plot_p_perm(perm_df)
perm_df %>%
arrange(desc(p_val))
For each of the 36 outcome variables, two mixed models (a base and an extended model) were fit and the extended model of interest was compared against the null/base model using a likelihood ratio test. The resulting p-values are reported. ICCs based on the extented model are also reported.
Of note, two of the models estimated a 0 variance for the site random effect. P-values and ICCs are not provided for these two models.
# Full models
run_full_models <- function(data){
models <- outcomes %>%
lapply(function(x) lmer(reformulate("(1|subject) + (1|site)", x), data=data)) %>%
setNames(outcomes)
}
# Null models
run_null_models <- function(data){
models_0 <- outcomes %>%
lapply(function(x) lmer(reformulate("(1|subject)", x), data=data)) %>%
setNames(outcomes)
}
#' Test variance term
#' @param mod full model
#' @param mod_0 null model
test_variance_term <- function(mod, mod_0){
A = lme4:::anova.merMod(mod, mod_0)
lrt_stat = as.numeric(-2*(logLik(mod_0) - logLik(mod)))
#if(near(lrt_stat, A$Chisq[2])) return(NA)
# Naive p-value given the chi_sq(df = 1) null distribution
pval_naive = pchisq(q = lrt_stat, df = 1, lower.tail = FALSE)
# p-value given a 50:50 mixture of chisq(1) and chisq(2),
# according to Self & Liang 1987
pval_SL = pchisq(q = lrt_stat, df = 1, lower.tail = FALSE)/2 +
pchisq(q = lrt_stat, df = 2, lower.tail = FALSE)/2
return(data.frame(pval_naive = pval_naive,
pval_SL = pval_SL))
}
# Which random effect variance is 0
get_zero_var_model_names <- function(models){
get_var <- function(model) summary(model)$varcor$site[1,1]
zero_var_models <- models %>%
purrr::map_dbl(get_var) %>%
as.data.frame() %>%
setNames('variance') %>%
rownames_to_column(var = 'outcome') %>%
filter(near(variance, 0, .Machine$double.eps)) %>%
pull(outcome)
return(zero_var_models)
}
models <- run_full_models(df)
models_0 <- run_null_models(df)
zero_var_model_names <- get_zero_var_model_names(models)
for (name in zero_var_model_names){
cat('###', name, '\n\n')
print(summary(models[[name]]))
cat("\n\n")
}
## ### MIMOSA_LESION_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -353.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5395 -0.3266 0.0607 0.3375 4.0443
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0013879 0.03725
## site (Intercept) 0.0000000 0.00000
## Residual 0.0004207 0.02051
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.32571 0.01147 28.39
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
make_SL_df <- function(models, models_0){
purrr::map2(models, models_0, test_variance_term) %>%
setNames(outcomes) %>%
bind_rows(.id = 'outcome') %>%
#filter(!(outcome %in% zero_var_model_names)) %>%
select(-pval_naive) %>%
rename(pval = pval_SL) %>%
mutate(pval_fdr = p.adjust(pval, method = 'fdr')) %>%
arrange(pval)
}
p_df <- make_SL_df(models, models_0)
plot_p_SL <- function(p_df){
p_df %>%
mutate(outcome = fct_reorder(outcome, pval_fdr, .desc = TRUE)) %>%
ggplot(aes(x=outcome, y=pval_fdr)) +
geom_point(size = 2) +
geom_hline(yintercept = 0.05, color = 'red') +
coord_flip() +
ggtitle("P-values of Likelihood Ratio Test (Self & Liang 1987)") +
ylim(0,1)
}
plot_p_SL(p_df)
options(scipen = 999)
p_df
ICC for the 35 remaining models is shown
options(scipen=999)
#' Get ICC for a single model
get_icc <- function(df){
outcomes %>%
map(~lmer(reformulate("(1|subject) + (1|site)", .x), data=df)) %>%
purrr::map(~ performance::icc(.x, by_group = TRUE) %>% as.data.frame) %>%
setNames(outcomes) %>%
bind_rows(.id = "outcome") %>%
select_if(~ !all(is.na(.))) %>% # eliminate column of NAs
na.omit() %>% # eliminate rows with NA
pivot_wider(names_from = 'Group', values_from = 'ICC', names_prefix = 'ICC_') %>%
relocate(ICC_site, .after = 'outcome') %>%
arrange(desc(ICC_site))
}
# extract observed icc_df
(icc_df <- get_icc(df))
## Make violin plots
## AND make line
plot_ICC_site <- function(icc_df, boot_df){
icc_df %>%
left_join(boot_df %>%
select(-ICC_site, -ICC_subject), by = 'outcome') %>%
mutate(outcome = fct_reorder(outcome, ICC_site)) %>%
ggplot(aes(x=outcome, y=ICC_site)) +
geom_errorbar(aes(ymin=lwr_site, ymax=uppr_site)) +
geom_point() +
coord_flip() +
ggtitle("ICC of Site Random Effect") +
ylim(0,1)
}
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)
plot_ICC_subject <- function(icc_df, boot_df){
icc_df %>%
left_join(boot_df %>%
select(-ICC_site, -ICC_subject), by = 'outcome') %>%
mutate(outcome = fct_reorder(outcome, ICC_subject)) %>%
ggplot(aes(x=outcome, y=ICC_subject)) +
geom_errorbar(aes(ymin=lwr_subject, ymax=uppr_subject)) +
geom_point() +
coord_flip() +
ggtitle("ICC of Subject Random Effect") +
ylim(0,1)
}
plot_ICC_subject(icc_df, boot_df)
After FDR-adjustment, 34 of the 36 variables show site effects according to the mixed models.
Here, we report results from the mixed models carried out only on the non-Hopkins data.
models <- df %>%
filter(site != 'Hopkins') %>%
run_full_models()
models_0 <- df %>%
filter(site != 'Hopkins') %>%
run_null_models()
zero_var_model_names <- get_zero_var_model_names(models)
for (name in zero_var_model_names){
cat('###', name, '\n\n')
print(summary(models[[name]]))
cat("\n\n")
}
## ### ATROPOS_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1067.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.55698 -0.51245 0.08047 0.54913 2.26173
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000002327184936992856181 0.0000152551136901
## site (Intercept) 0.0000000000000000000000008819 0.0000000000009391
## Residual 0.0000000000543621244491441043 0.0000073730675061
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000757099 0.000004715 160.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1067.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.69258 -0.49167 0.02313 0.54945 2.34842
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000269168850518 0.000016406366
## site (Intercept) 0.000000000000000007242 0.000000002691
## Residual 0.000000000052691491108 0.000007258890
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000695292 0.000005051 137.7
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00439263 (tol = 0.002, component 1)
p_df <- make_SL_df(models, models_0)
plot_p_SL(p_df)
options(scipen = 999)
p_df %>%
arrange(pval_fdr)
ICC for the 35 remaining models is shown below
# extract icc_df
icc_df <- df %>%
filter(site != 'Hopkins') %>%
get_icc()
icc_df
options(scipen=999)
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)
plot_ICC_subject(icc_df, boot_df)
After FDR-adjustment, 14 of the 36 variables show site effects according to the mixed models run on non-Hopkins data.
Here, we report results from the mixed models carried out only on the non-Hopkins data.
models <- harmonized_df %>%
run_full_models()
models_0 <- harmonized_df %>%
run_null_models()
zero_var_model_names <- get_zero_var_model_names(models)
for (name in zero_var_model_names){
cat('###', name, '\n\n')
print(summary(models[[name]]))
cat("\n\n")
}
## ### ATROPOS_GM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1470.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9497 -0.6559 0.1219 0.5873 1.9042
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000018824 0.000013720
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000006596 0.000008122
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000875530 0.000004252 205.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1510.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5691 -0.4094 0.0438 0.4148 2.2466
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000011324 0.000010641
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000003757 0.000006129
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000891003 0.000003293 270.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1610.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7396 -0.6439 0.1238 0.6573 1.9164
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000012776 0.000011303
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000005367 0.000007326
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000865351 0.000003509 246.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1635.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5691 -0.4683 0.0466 0.4440 1.9366
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000011448 0.000010700
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000003799 0.000006163
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000898376 0.000003302 272.1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1493.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.89528 -0.56640 0.04294 0.66129 1.63512
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000564024297156820700447 0.000023749195716
## site (Intercept) 0.000000000000000000000000000121 0.000000000000011
## Residual 0.000000000237605972500836866663 0.000015414472826
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000881333 0.000007374 119.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1581.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2746 -0.6955 0.1126 0.7008 1.6449
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000043022 0.00002074
## site (Intercept) 0.00000000000000 0.00000000
## Residual 0.00000000006889 0.00000830
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000918269 0.000006325 145.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1554.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1783 -0.5615 -0.0164 0.6371 1.9517
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000003024 0.00001739
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000001069 0.00001034
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000972874 0.000005375 181
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1499.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.91188 -0.58387 0.09154 0.69656 1.55554
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000005985 0.00002446
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000002165 0.00001472
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000873548 0.000007565 115.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1307.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2819 -0.6582 0.0197 0.4861 2.4428
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000001033819770882159319180 0.0001016769281048
## site (Intercept) 0.00000000000000000000000007073 0.0000000000002659
## Residual 0.00000000230136384853161502135 0.0000479725322297
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00124147 0.00003114 39.86
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -615.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5514 -0.5452 -0.0296 0.4591 4.5337
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000039622 0.006295
## site (Intercept) 0.000000000 0.000000
## Residual 0.000007389 0.002718
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.160736 0.001926 83.45
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -534.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1642 -0.5512 0.0460 0.5211 2.1663
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00043293 0.02081
## site (Intercept) 0.00000000 0.00000
## Residual 0.00001823 0.00427
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.387185 0.006295 61.51
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -642.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1355 -0.4601 -0.0775 0.5206 4.2706
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00005969 0.007726
## site (Intercept) 0.00000000 0.000000
## Residual 0.00001007 0.003174
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.173234 0.002357 73.48
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -572.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.99637 -0.54746 0.06809 0.49398 2.41775
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00056946 0.023863
## site (Intercept) 0.00000000 0.000000
## Residual 0.00002005 0.004478
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.403270 0.007213 55.91
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -526.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92394 -0.68760 0.09937 0.63262 2.03870
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00036149 0.019013
## site (Intercept) 0.00000000 0.000000
## Residual 0.00004206 0.006485
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.30928 0.00578 53.51
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -690.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67793 -0.65587 0.09039 0.70079 1.91229
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000042752 0.006539
## site (Intercept) 0.000000000 0.000000
## Residual 0.000005261 0.002294
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.149733 0.001989 75.29
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -473.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78535 -0.67347 -0.08739 0.59236 3.04956
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00025223 0.01588
## site (Intercept) 0.00000000 0.00000
## Residual 0.00009467 0.00973
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.319816 0.004916 65.06
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -350.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2050 -0.3210 0.0810 0.3787 3.7290
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.001507 0.03882
## site (Intercept) 0.000000 0.00000
## Residual 0.000434 0.02083
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.32570 0.01195 27.27
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1483.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2608 -0.5519 0.0936 0.4934 2.1900
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000022050 0.000014849
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000005204 0.000007214
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000753768 0.000004561 165.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1596
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9838 -0.4482 0.1488 0.5683 1.6362
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000241102 0.000015527
## site (Intercept) 0.000000000000000 0.000000000
## Residual 0.000000000008663 0.000002943
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000615973 0.000004695 131.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1632.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3337 -0.5809 0.0119 0.4674 2.4205
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000017616 0.000013272
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000003744 0.000006118
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000736281 0.000004062 181.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1715.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0897 -0.5025 0.1614 0.7020 1.6303
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000027296 0.000016522
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000001059 0.000003254
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000611652 0.000004995 122.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1570.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.92049 -0.48499 -0.00084 0.59701 2.01372
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000020565 0.000014341
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000008974 0.000009473
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000663841 0.000004457 148.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1605.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00814 -0.67017 0.06203 0.68540 1.95853
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000049202 0.000022182
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000004805 0.000006932
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000802454 0.000006735 119.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1652.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.65999 -0.63934 0.01372 0.67020 1.81063
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000027771 0.000016665
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000002641 0.000005139
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000649447 0.000005059 128.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1569.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0016 -0.5577 0.0336 0.6803 1.8869
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000030174 0.000017371
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000008641 0.000009296
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000651171 0.000005344 121.8
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1314.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9024 -0.5410 -0.0192 0.4109 3.2522
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000005411 0.00007356
## site (Intercept) 0.000000000000 0.00000000
## Residual 0.000000002272 0.00004766
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00091150 0.00002284 39.91
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1484.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3266 -0.5036 0.0753 0.5078 2.2200
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000025272 0.000015897
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000005024 0.000007088
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000692596 0.000004869 142.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1602.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39757 -0.52379 0.01794 0.64805 2.18877
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000407377 0.000020184
## site (Intercept) 0.000000000000000 0.000000000
## Residual 0.000000000007204 0.000002684
## Number of obs: 74, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000478127 0.000006094 78.45
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1634.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5475 -0.4727 -0.0315 0.5050 2.5274
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000022075 0.000014858
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000003514 0.000005928
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000671485 0.000004531 148.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1712.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.98135 -0.72897 0.07567 0.74622 2.14054
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000049321 0.000022208
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000001013 0.000003182
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000467961 0.000006706 69.78
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1604.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.71875 -0.49645 -0.00331 0.50080 2.37664
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000022539 0.000015013
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000005451 0.000007383
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000554630 0.000004605 120.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1613.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9257 -0.6679 0.0833 0.6353 2.0761
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000053911 0.000023219
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000004216 0.000006493
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00074452 0.00000704 105.8
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1704.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9020 -0.6091 -0.1556 0.6481 2.1291
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000462825961584492 0.00002151339029
## site (Intercept) 0.000000000000000000007618 0.00000000008728
## Residual 0.000000000011480548228789 0.00000338829577
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000487319 0.000006498 74.99
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1580.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.55747 -0.62126 0.09323 0.66516 1.91966
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000028167 0.000016783
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000007436 0.000008623
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000539506 0.000005155 104.7
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1309.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0764 -0.5268 -0.0541 0.3975 3.4756
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000005008 0.00007077
## site (Intercept) 0.000000000000 0.00000000
## Residual 0.000000002472 0.00004972
## Number of obs: 80, groups: subject, 11; site, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00074756 0.00002208 33.86
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
p_df <- make_SL_df(models, models_0)
plot_p_SL(p_df)
options(scipen = 999)
p_df %>%
arrange(pval_fdr)
ICC for the 35 remaining models is shown below
# extract icc_df
icc_df <- harmonized_df %>%
get_icc()
icc_df
options(scipen=999)
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_harmonized_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)
plot_ICC_subject(icc_df, boot_df)
After FDR-adjustment, 14 of the 36 variables show site effects according to the mixed models run on non-Hopkins data.
Here, we report results from the mixed models carried out only on the non-Hopkins data.
models <- harmonized_no_hopkins_df %>%
run_full_models()
models_0 <- harmonized_no_hopkins_df %>%
run_null_models()
zero_var_model_names <- get_zero_var_model_names(models)
for (name in zero_var_model_names){
cat('###', name, '\n\n')
print(summary(models[[name]]))
cat("\n\n")
}
## ### ATROPOS_GM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1057.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59679 -0.56995 -0.00606 0.47152 2.21326
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000019723 0.000014044
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000007148 0.000008455
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000880281 0.000004398 200.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1123.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.56177 -0.61044 0.04468 0.46739 1.89250
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000011658 0.000010797
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000001751 0.000004185
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000898208 0.000003308 271.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1197.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59468 -0.54732 -0.07949 0.49328 2.18728
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000012838 0.000011331
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000005379 0.000007334
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000869450 0.000003549 245
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1250.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67798 -0.51885 -0.02254 0.49969 1.87665
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000001268 0.000011260
## site (Intercept) 0.0000000000000 0.000000000
## Residual 0.0000000000184 0.000004289
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000904548 0.000003441 262.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1108.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70139 -0.41705 -0.00096 0.41623 1.73233
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000007067 0.00002658
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000002326 0.00001525
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000881549 0.000008261 106.7
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1170.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0924 -0.5592 -0.1059 0.6010 1.7205
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000042429 0.000020598
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000007354 0.000008575
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000924270 0.000006312 146.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1165.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.48121 -0.50891 -0.05478 0.59726 1.96538
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000036606 0.000019133
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000008413 0.000009172
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000977541 0.000005893 165.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1111.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.74833 -0.43997 -0.02685 0.49318 1.72483
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000007461 0.00002732
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000002162 0.00001470
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000873240 0.000008458 103.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_AD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_AD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1005.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25114 -0.66721 -0.04421 0.52226 2.38080
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000012896 0.00011356
## site (Intercept) 0.000000000000 0.00000000
## Residual 0.000000001061 0.00003257
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00125981 0.00003451 36.51
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -434.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9565 -0.5071 -0.0176 0.4376 3.9716
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000047638 0.006902
## site (Intercept) 0.000000000 0.000000
## Residual 0.000007877 0.002807
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.161855 0.002118 76.41
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -397.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06094 -0.58666 0.06531 0.60149 2.30536
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00046197 0.021493
## site (Intercept) 0.00000000 0.000000
## Residual 0.00001105 0.003324
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.390175 0.006497 60.05
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -465.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8100 -0.4464 -0.0713 0.4812 3.7837
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00006570 0.008105
## site (Intercept) 0.00000000 0.000000
## Residual 0.00001136 0.003371
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.174061 0.002484 70.09
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -420.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.97762 -0.50314 -0.00207 0.51873 2.26368
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00058732 0.02423
## site (Intercept) 0.00000000 0.00000
## Residual 0.00001832 0.00428
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.404517 0.007329 55.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -385.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7845 -0.6547 0.1306 0.6505 2.0661
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00041177 0.020292
## site (Intercept) 0.00000000 0.000000
## Residual 0.00003988 0.006315
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.309028 0.006174 50.05
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -450.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6524 -0.6892 0.1759 0.5564 2.0740
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00048879976665297715250 0.0221088164915
## site (Intercept) 0.00000000000000000001728 0.0000000001314
## Residual 0.00001037093894721858281 0.0032203942223
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.428463 0.006679 64.15
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -341.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6245 -0.6431 -0.1295 0.6183 2.5472
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0003010 0.01735
## site (Intercept) 0.0000000 0.00000
## Residual 0.0001039 0.01020
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.321027 0.005399 59.46
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_FA
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -285.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1851 -0.5027 0.1120 0.4741 4.0323
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0018111 0.04256
## site (Intercept) 0.0000000 0.00000
## Residual 0.0002298 0.01516
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.32520 0.01298 25.05
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1061.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.74244 -0.50704 -0.00033 0.40066 2.44188
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000023279 0.000015258
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000006331 0.000007957
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000757206 0.000004734 159.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1175.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.34089 -0.42407 0.00377 0.57602 2.29103
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000242658 0.00001558
## site (Intercept) 0.000000000000000 0.00000000
## Residual 0.000000000004408 0.00000210
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000619108 0.000004706 131.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1204.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.02940 -0.44962 -0.04102 0.45851 2.65693
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000017208 0.000013118
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000004435 0.000006659
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00073933 0.00000405 182.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1302.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2217 -0.6082 0.2078 0.6628 1.9767
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000273473 0.000016537
## site (Intercept) 0.000000000000000 0.000000000
## Residual 0.000000000005432 0.000002331
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000615070 0.000004995 123.1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1163.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70260 -0.41882 -0.09437 0.58217 1.94969
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000023837 0.000015439
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000009383 0.000009687
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000664123 0.000004825 137.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1182.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.17470 -0.61022 0.01029 0.64407 1.97465
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000048376 0.00002199
## site (Intercept) 0.00000000000000 0.00000000
## Residual 0.00000000005639 0.00000751
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000806894 0.000006704 120.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1227
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.35390 -0.53871 0.00628 0.50636 1.74888
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000028359 0.000016840
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000002528 0.000005028
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00065111 0.00000512 127.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1155.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.80432 -0.45017 0.02758 0.58016 2.04464
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000003319 0.00001822
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000001045 0.00001022
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000650186 0.000005654 115
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_MD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_MD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1022.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.04787 -0.60690 0.02819 0.50027 2.43177
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000073893 0.00008596
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000008431 0.00002904
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.0009252 0.0000262 35.32
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1059.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.74161 -0.42889 0.02642 0.47831 2.42025
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000002688303057123714644 0.0000163960454291
## site (Intercept) 0.0000000000000000000000002266 0.0000000000004761
## Residual 0.0000000000637590763813453111 0.0000079849280762
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000695413 0.000005069 137.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### ATROPOS_WM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1166
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.30250 -0.55062 -0.08627 0.62789 2.13673
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000421592 0.000020533
## site (Intercept) 0.000000000000000 0.000000000
## Residual 0.000000000004868 0.000002206
## Number of obs: 54, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000479304 0.000006199 77.32
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1201.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1286 -0.4092 -0.0490 0.4789 2.7077
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000216717387668859 0.00001472132425
## site (Intercept) 0.000000000000000000005284 0.00000000007269
## Residual 0.000000000045005642021951 0.00000670862445
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000674042 0.000004525 149
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FAST_WM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1274.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9017 -0.5845 0.0928 0.6056 2.0575
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.000000000501142 0.000022386
## site (Intercept) 0.000000000000000 0.000000000
## Residual 0.000000000008631 0.000002938
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000470076 0.000006761 69.53
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### FIRST_THALAMUS_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1184.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.44872 -0.64361 -0.08354 0.52868 1.98840
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000002388 0.000015452
## site (Intercept) 0.0000000000000 0.000000000
## Residual 0.0000000000622 0.000007887
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000555019 0.000004772 116.3
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_GM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1184.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.10833 -0.54398 0.01153 0.63777 2.01162
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000005320 0.000023065
## site (Intercept) 0.0000000000000 0.000000000
## Residual 0.0000000000531 0.000007287
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00074820 0.00000702 106.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_WM_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1251.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75003 -0.46773 0.00206 0.46903 1.90570
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.00000000046178 0.000021489
## site (Intercept) 0.00000000000000 0.000000000
## Residual 0.00000000001385 0.000003721
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000487561 0.000006498 75.04
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### JLF_THALAMUS_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1156.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3887 -0.6660 0.1540 0.6069 1.8568
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000002898 0.00001702
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000001041 0.00001020
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.000538321 0.000005304 101.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
##
##
##
## ### MIMOSA_LESION_RD
##
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_RD ~ (1 | subject) + (1 | site)
## Data: data
##
## REML criterion at convergence: -1021.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2834 -0.7032 0.0896 0.6080 2.2461
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.0000000069865 0.00008359
## site (Intercept) 0.0000000000000 0.00000000
## Residual 0.0000000008635 0.00002939
## Number of obs: 60, groups: subject, 11; site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.00075895 0.00002549 29.77
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
p_df <- make_SL_df(models, models_0)
plot_p_SL(p_df)
options(scipen = 999)
p_df %>%
arrange(pval_fdr)
ICC for the 35 remaining models is shown below
# extract icc_df
icc_df <- harmonized_no_hopkins_df %>%
get_icc()
icc_df
options(scipen=999)
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_harmonized_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)
plot_ICC_subject(icc_df, boot_df)
After FDR-adjustment, 14 of the 36 variables show site effects according to the mixed models run on non-Hopkins data.